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Abstract. We applied the simulated tempering and magnetizing (STM) method 
to the two-dimensional three-state Potts model in an external magnetic field 
in order to perform further investigations of the STM's applicability. The 
temperature as well as the external field are treated as dynamical variables 
updated during the STM simulations. After we obtained adequate information 
for several lattice sizes L (up to 160 X 160), we also performed a number of 
conventional canonical simulations of large lattices, especially in order to illustrate 
the crossover behavior of the Potts model in external field with increasing L. The 
temperature and external field for larger lattice size simulations were chosen by 
extrapolation of the detail information obtained by STM. We carefully analyzed 
the crossover scaling at the phase transitions with respect to the lattice size as well 
as the temperature and external field. The crossover behavior is clearly observed 
in the simulations in agreement with theoretical predictions. 
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1. Introduction 

Monte Carlo (MC) and molecular dynamics (MD) methods have been demonstrated 
in many applications to be indispensable tools to study the statistical properties 
of various systems in equilibrium. The quasi-ergodicity problem, however, where 
the system gets trapped in states of local energy minima, has often posed great 
difficulties. In order to overcome this difficulty, generalized-ensemble algorithms have 
been developed and applied to many problems including spin models and biomolecular 
systems (for reviews, see, e.g., Refs. [I, 2, 3, 4]). 

Well-known examples of generalized-ensemble algorithms are the multicanonical 
algorithm (MUCA) [5, 6], simulated tempering (ST) [7, 8], and the replica-exchange 
method (REM) [9, 10] (also referred to as parallel tempering). Closely related to 
MUCA arc the Wang-Landau method [II, 12] and metadynamics [13]. REM is 
implicitly a special case of the general method described in the earlier work of Ref. 
[14], as detailed later in Ref. [15]. 

Based on the recent multi-dimensional generalization of generalized-ensemble 
algorithms [16, 17, 18], the "Simulated Tempering and Magnetizing" (STM) method 
has been proposed and developed [19, 20]. In Refs. [19, 20] two of us studied the 
classical Ising model, introducing the external (magnetic) held as a second dynamical 
variable besides the temperature and showed improvements over the conventional 
"one-dimensional" simulated tempering schemes, such as better sampling efficiency 
and potential applicability to a first-order phase transition, which cannot be dealt 
with by one-dimensional ST. 

In the present work, we further investigate the STM method, applying it to the 
two-dimensional three-state Potts model in external magnetic field [21, 22] . This model 
has several interesting applications in condensed matter physics [22] and its three- 
dimensional counterpart serves as an effective model for quantum chromodynamics 
[23, 24, 25, 26]. We see the STM scheme working in this more complicated system 
as well. We also look into crossover behaviors according to lattice size L as well 
as temperature T and external field h. We observe that the STM method lets us 
investigate a wide area of sampling space. 

The rest of this article is organized as follows. In section 2 we review the STM 
method and give the details of our simulations. In section 3 we present the results. 
After checking the two-dimensional random walks, we compare ST and STM, and 
calculate various thermodynamic quantities at many sets of parameter values, in 
combination with reweighting techniques. We then use this extensive data set to study 
the crossover scaling behavior at the phase transitions with respect to the lattice size 
as well as the temperature and external field. In section 4 we conclude this article. 
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2. Model and Methods 

2.1. Model 

We study the two-dimensional three-state standard Potts model in external field with 
energy 

H = E-hM, (1) 

B = -£<w,. ( 2 ) 

(id) 

N 

M = ^V 4 > (3) 

i=l 

where N = L 2 denotes the total number of spins, 5 is the Kronecker delta function, cr, 
a spin at the i-th site, and ft. the external field. The spin Oi takes on one of the three 
values 0, 1, or 2. The sum in (2) goes over all nearest-neighbor pairs, where the spins 
are arranged on a square L x L lattice with periodic boundary conditions. Data were 
obtained in STM for lattice sizes ranging from 2 x 2 to 160 x 160 and additionally in 
conventional canonical simulations on 320 x 320 and 640 x 640 lattices. 

We recall that the three-state (standard) Potts model is equivalent to the three- 
state planar Potts or Z 3 model. We first introduce a spin 

/cos ^crA 

— U*J- (4 » 



The zero-field energy term is then given by 



and the magnetization reads 



-^planar — ^ ^ s i ' s j > (5) 



N 

M = $>. (6) 



In an external field chosen along the x-direction, 

h = he x , (7) 
the product of external field and magnetization becomes 

h-M = hM^\ (8) 

where stands for the x-component of M, which is given by ^ cos ^Cj. Because 
the argument of the cosine and sine in (4) can only take the values of 0, 27r/3, and 
47r/3, we have 

E = (£ pW - JV) x H , ( 9 ) 
M = (m& + y ^ x I . (10) 
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Thus, we arrive at 

2 2 , v N 

E - hM = -£ planar - -WW - —(2 - h) . (11) 

Because the last term is an unimportant constant, the standard Potts model 
is equivalent to the planar Potts model with a | -smaller coupling constant and 
magnetization normalized by |: 

H = E-hM, (12) 
£ = -H^cos%, (13) 

2 N 

M=-Vcos^, (14) 
6 i=i 

where 0y and #j are defined by 6j — 0j and ^-cr.j, respectively. 

Because M' x ' is the projection of M onto the x-axis, M^/N equals 1, -1/2, 
— 1/2, and when the system is in the O-direction ordered, the 1-direction ordered, 
the 2-direction ordered, and the disordered state, respectively. This corresponds to 
M/N — when the system is in 1- or 2-direction ordered states and M/N = 1/3 in 
the disordered state. 

As one can see from (1) and (3) [or, (12) and (14)], spin direction is favored by a 
positive external field (h > 0). Accordingly, a negative external field (h < 0) disfavors 
spin direction 0. Thus, the system is expected to behave like a two-dimensional Ising 
model in the presence of a negative external field. In fact, in the limit h — oo, 
the three-state Potts model is equivalent to the Ising model in zero external field, 
because the unfavored states disappear in the partition function calculations. Figure 
1 illustrates the schematic picture of this relation. 




Figure 1. Schematic description of the behavior of spins according to the external 
field, (a) Spin is favored at h > 0, i.e., spin 1 and spin 2 are disfavored, (b) All 
the three states are equivalent at h = 0. (c) Spin is disfavored at h < 0. 



For the easier availability of reference, we summarize the critical exponents for 
the two-dimensional Ising and three-state Potts model in table 1 [22]. 

2.2. Simulation methods 

In this section we briefly review the STM method [19, 20]. While the conventional 
ST method [7, 8] considers the temperature to be the dynamical variable, the STM 
method employs the external field as a second dynamical variable besides temperature. 
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Table 1. Critical exponents for the two-dimensional Ising and three-state Potts 



models (yt = 


Vh = 


03 + 


7)/")- 






Model yt 


Vh 




7 




V 


Ising 1 


15/8 


1/8 


7/4 


15 


1 


Potts 6/5 


28/15 


1/9 


13/9 


14 


5/6 



This algorithm is based on the multi-dimensional extension of generalized-ensemble 
algorithms [16, 17, 18]. In other words, we consider 

e -(E-hM)/T+a(T,h) ( 15 ) 

as a joint probability for (x, T, h) (eX® {Ti, T 2 , . . . , Tjv T } <8> {hi, hi, ■ ■ ■ fiN h }) 7 where 
a(T,h), x, and X are a parameter, the (microscopic) state, and the sampling space, 
respectively. Here and hereafter, we set Boltzmann's constant to unity. Note that the 
temperature and external field are discretized into Nt and Nh values, respectively. 

We will find the candidate for a(Ti,hj) by looking into the probability of 
occupying each set of parameter values. It is given by 

f E(x)-h, M(x) , ,_ , . 

P{Tu hi) = (16) 



fe=l 1=1 



EE/ da;e " 

e -/(T i ,/i j )+o(T i ,/i j ) 

y^y" e -/(T*,hi)+<»(r*,/n) 
fc=i ;=i 

-/(T^h^+oCTi,^) 



(17) 



where 

e -m, hj) = f dxe -(E- hj M)/T K (19) 

Thus, the dimensionless free energy /(Ti, /ij) is a proper choice for a(Ti, hj) in order to 
acquire a uniform distribution of the number of samples according to T and h. These 
values can be estimated by a number of methods. For example, one can obtain such 
values from preliminary simulations and reweighting techniques. 

Any thermal average (A) Tih at given Tj (g {Ti, T2, . . . , Tn t }) and hj (s 
{fti, hi, • • • , hN h }) can be obtained by the conditional expectation: (A) T . h . = 
(A\Ti,hj) STM . Namely, we have 



N Ti,h } 

( A h, hj = J2 A k hj . (20) 

*' 3 k=l 

where Nx^hj is the total number of samples obtained at T, and hj, and A\, h . 
represents the fc-th sample at Tj and hj. 
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Table 2. 


Parameter values for the STM simulations. 








Lattice size L 


5 


10 


20 


40 


80 


160 


Total MC sweeps 


220915510 


216277330 


164000000 


164000000 


43605000 


525000000 


Parameter-updating period 


1 


1 


1 


1 


1 


1 


Ti , . . . , Tj\r T 


0.1, ... ,3.5 


0.1, ...,3.5 


0.3,..., 3. 2 


0.3, ... ,3.2 


0.5, ...,2.0 


0.5, ... ,2.0 


hi, ... , h Nh 


-1.5,..., 1.5 


-1.5,..., 1.5 


-1.5, ...,1.5 


-1.5, ... ,1.5 


-1.5, ...,1.5 


-1.0, . . . ,0.025 


N T 


20 


20 


40 


40 


76 


95 


N h 


21 


21 


41 


41 


5 1 


30 


A^data " 


10 


10 


10 


10 


10 


10 



The data were stored every -/V<i ata MC sweeps. 



The temperature T or external field h can be updated similarly to the spin <Tj, 
because they are considered as dynamical variables. The Metropolis criterion for 
updating T and h is given by 

w(T it h d -+ , hj> ) = min (l, ^r^y) (21) 

= min(l,exp(- (±- - E + Qj- - M + aft,, h r ) - a(T u fc,-))) . 

(22) 

Once an initial state is prepared, the STM simulations can be performed by 
repeating the following two steps: 1. We perform a conventional canonical simulation 
at fixed T and hj for a certain number of MC sweeps. 2. We update the temperature 
and/or external field by (22) with a(T, h) = f(T, h). 

In our implementation, every fixed number of MC sweeps, either T or h was 
updated (the choice between T and h was made at random) by (22) to a neighboring 
value (the choice between two possible neighbors was also made at random). Here, one 
MC sweep consists of L x L single spin updates. The number of MC sweeps performed 
between parameter updates is here referred to as the parameter-updating period. 

We remark that, as spins can be updated by a number of algorithms, other 
schemes of updating the parameters can be employed [27]. There also exists a 
temperature updating scheme for ST by the Langcvin algorithm [28] . 

Table 2 summarizes the conditions of the present STM simulations. According to 
the previous studies [20, 29, 30], we updated the parameters frequently. That is, we 
employed very small parameter-updating periods. 

In addition, we also performed conventional canonical simulations. Table 3 lists 
their details. The temperature was chosen by extrapolations of the STM results. We 
estimated the proper temperature by fitting the STM results to T max — T c oc L~ x l v , 
where T max is the temperature at which the observables take their maxima. The Greek 
letter v denotes the critical exponent of the correlation length. For vanishing external 
field we fitted the data to the Potts case (y = 5/6) and in (negative) external field to 
the Ising case (v = 1), respectively. 

As for spin updates, we employed the single-spin update algorithm; we updated 
spins one by one with the heatbath algorithm. As for quasi-random-number generator, 
we used the Mersenne Twister [31]. 
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Table 3. Parameter values for the regular canonical simulations. 



Lattice size L 320 640 



Total MC sweeps a 4000000 4000000 

(T, h) b (0.995518, 0), (0.995995, 0), (0.996490, 0) (0.994985, 0), (0.995209, 0), (0.995512, 0) 

(1.0752077, -0.5), (1.077044, -0.5), (1.078302, -0.5) (1.075672, -0.5), (1.076181, -0.5), (1.076936, -0.5) 

(1.101447, -1.0), (1.102321, -1.0), (1.103380, -1.0) (1.100386, -1.0), (1.101043, -1.0), (1.101811, -1.0) 

^data c 10 20 

a The number performed for each set of temperature and external field. 
b Three temperature values for each external field. 
c The data were stored every Af^ata MC sweeps. 



2.3. Free energy calculations 

The simulated tempering parameters, or free energy, in (15) and (19) can be simply 
obtained by reweighting techniques applied to the results of preliminary simulation 
runs [16, 17, 18, 32]. We used two reweighting methods for the free energy calculations. 
One method is the multiple-histogram reweighting method, or Weighted Histogram 
Analysis Method (WHAM) [33, 34, 35] and the other is the Multistate Bennett 
Acceptance Ratio estimator (MBAR) [36], which is based on WHAM. 

The equations of the WHAM algorithm applied to the system are as follows. 
For details, the reader is referred to Refs. [17, 34, 35]. The density of states (DOS) 
n(E, M) and free energy values f(T i7 hj) can be obtained from 

N T N h 

n (E, M) = -j^~jf h ^ , (23) 

E E N ^ ex P(/( T - h i) -i E ~ hjM)/^) 
»=i j=i 

f(Ti,hj) - - In E M ) exp(-(B - hjM)/TJ , (24) 

E,M 

where n^,^ M) is the histogram of E and M at Tj and hj, and Nx^hj is the 
total number of samples obtained at T, and hj. By solving these two equations self- 
consistently by iteration, we can obtain n(E, M) and /(Tj, hj). The obtained n(E, M) 
allows one to calculate any thermal average at arbitrary temperature and external field 
values. Note that f(Ti,hj) is determined up to a constant, which sets the zero point 
of free energy. Accordingly, n(E, M) is determined up to a normalization constant. 

The MBAR is based on the following equations. Namely, by combing (23) and 
(24) , the free energy can be written as 

"S tsts "rt-'^-Mwr.) , (25) 



n=l 

fe=l 1=1 



E E N n,hi exp(/(T fc) hi) - (E n - h t M n )/T k ) 



where N, Nx k ,h u E n , and M n is the total number of data, the number of samples 
associated with and hi, energy of the n-th data, and magnetization of the n-th 
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data, respectively. This equation should be solved self-consistently for f(T it hj). Note 
that, as in WHAM, f(Ti,hj) is determined up to a constant. 

We repeat the preliminary STM simulations and free energy calculations until we 
finally obtain sufficiently accurate free energy values which let the system perform a 
random walk in the temperature and external field space during the STM simulation. 
We then perform a final production run. 

Note that these two reweighting methods enable us to obtain not only 
dimensionless free energy values but also physical values at any temperature and at 
any external field. Such averages are given by 

N 

( A )T,h = 52W na A{x n ), (26) 

n=l 

w 1 exp(-(E n -hM n )/T) 

Wna ~ ~^r h ' [Z1 > 



J2 E N n,hi ex P (/(T fe , hi) - {E n - hiM n )/T k ) 



fe=i i=i 

N 



ln S_ST cM-(E n -hM n )/T) 

[Ca) ~ 2^ Nt Nh • \ Zti ) 



n=l 

fe=l 1 = 1 



J2 N n,hi exp(/(T fe , hi) - (E n - hiM n )/T k ) 



For details, the reader is referred to Refs. [36, 37]. 

We also used two other methods of free energy calculations. One is given as 
follows. By substituting a(T,h) in (18) by the estimates for free energy f(T,h), we 
obtain 

P(T,ft)«e-^W"( T ' k ). (29) 

We can write 

f(T, h) = f(T, h) - In P(T, h) + const. (30) 

Here, P(T, h) can be obtained as a histogram at each set of parameter values in a 
preliminary STM simulation. Thus, this equation enables one to refine the free energy 
much more easily than the reweighting methods, because the method does not require 
any iterations. This method does not work well, however, when P(Tj, hj) is too small 
(or /(Tj, hj) is too far away from the true values) to obtain samples at (Tj, hj), while 
the reweighting techniques still work. 

The other method for the free energy calculations is a Wang-Landau-like scheme, 
where we subtract the free energy value being sampled by a fixed constant during 
preliminary simulations. To stand on the safe side, we did not use such data for 
reweighting techniques which, strictly speaking, requires equilibrium data as inputs. 
Note that this method also works with inaccurate free energy values. Thus, this 
method works even when the free energy estimates are far away from sufficiently 
accurate values. 

In the present work, we first used the reweighting methods and Wang-Landau-like 
scheme to obtain rough estimates of the free energy for the entire parameter space. We 
then used the combination of the reweighting methods and (30) for further refinements 
of the free energy. 
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We first examine whether the STM simulations were carried out properly or not. 
Figure 2 and figure 3 show the time series of temperature and external field, 
respectively, for L = 80. In both plots we see block structures reflecting the first- 
order phase transition line at h = in the Potts model (see figure 3) and the second- 
order phase transition at the effective Ising transition temperature T c (h) « 1.1346 
for negative external field (see figure 2). Within these blocks, the temperature and 
external field indeed realized random walks. 




0.25 0.5 0.75 1 

x10 6 MCsweep 



Figure 2. History of temperature T for the linear lattice size L = 80. 




0.25 0.5 0.75 

x10 6 MCsweep 



Figure 3. History of external field h for the linear lattice size L = 80. 

Figure 4 and figure 5 show the energy per spin and the magnetization per spin, 
respectively, as functions of MC sweeps. They also performed random walks. Note 
that there exist expected correlations between the temperature and energy (see figures 
2 and 4) and between the external field and magnetization (see figures 3 and 5). The 
same behavior was observed in simulations with other lattice sizes (data not shown). 
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0.25 0.5 0.75 1 
x10 6 MCsweep 



Figure 4. History of energy per spin, E/L 2 , for the linear lattice size L = 80. 



0.5 



0.25 




0.75 



x10 MCsweep 



Figure 5. History of the magnetization per spin, M/L 2 , for the linear lattice size 
L = 80. Recall that M/L 2 = 1 corresponds to the ordered state in 0-direction, 
M/L 2 = to the ordered states in 1- or 2-direction, and M/L 2 = 1/3 to the 
disordered state. 



Figure 6 shows the history of a differently defined magnetization given by 

3 



max 

j=0,l,2 



L 




Hereafter, we also use the following definition: 



L 2 



(31) 



(32) 



This quantity takes physically more intuitive values of 1 and when the system is in 
one of the three ordered phases and in the disordered phase, respectively. Here, we 
see a clear negative correlation between E/L 2 and M max /L 2 (see figures 4 and 6). 

To compare the results with ordinary ST simulations, we also performed a ST 
simulation with L = 40. The ST simulation was performed at the conditions similar 
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to those of STM; namely, the same total number of MC sweeps, same temperature 
distribution, and so on, except that here we set h = 0. 

With the data obtained, we performed the WHAM calculations to obtain the 
DOS. As shown in figure 7, the area sampled by STM is larger than that by ST. Thus, 
the STM method enables us to perform reweighting techniques in a wider range. We 
recall that M is zero in the 1- or 2-direction ordered phases. The disordered phase 
corresponds to M/L 2 = |. 



(a! x10 . 3 (b) x10 „ 3 




-2 -1.6 -1.2 -0.8 -2 -1.6 -1.2 -0.8 



E/L 2 E/L Z 

Figure 7. Calculated DOS obtained by WHAM with (a) ST and (b) STM data. 
The linear lattice size L was 80. 

We further closely look into the difference in the sampled areas between the 
two methods. Figure 8 illustrates how the sampled areas differ. The red, green, 
blue, and white regions correspond to the region that was sampled by the STM 
method exclusively, by both methods, by the ST method only, and by neither of 
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them, respectively. Thus, at first sight, it seems that there are some areas in which 
STM is not good and that ST is somehow more powerful than STM. 




E/L 2 



Figure 8. Difference in sampled points between ST and STM. The red, green, 
blue, and white regions correspond to the area sampled by only STM, by both of 
them, only by ST, and by neither of them, respectively (L = 40). 





Figure 9. Difference in sampled points between ST and STM. The red, green, 
blue, and white regions correspond to the area sampled by only STM, by both 
of them, only by ST, and by neither of them, respectively. Plots (a) and (b) are 
with different focus. 



Figure 9 zooms in on a region where blue is dominant (mainly sampled by 
ST). There are many pigments (in red and green) which both methods sampled and 
which even only STM sampled. This shows that because the ST method has more 
samples at a smaller number of parameter values, the part sampled is narrower but 
denser. However, the representative parts should be sampled properly by STM as 
well, although the sample density decreases. 
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To make it sure that the STM method also samples the relevant area sufficiently, 
we then performed reweighting analyses along h = with data obtained by ST 
and STM. Figure 10 and figure 11 show specific heat C/L 2 and susceptibility x/L 2 , 
respectively, as functions of T along the line h = 0. They are defined by 

(E 2 )-(E) 2 , 
C = K ' T2 K ' , (33) 

X = j, • (34) 

The red and green curves correspond to the data obtained by STM and by ST, 
respectively. The error bars were obtained by the jackknife method [38, 39, 40, 41]. 
We see no outstanding differences between the two methods. Thus, we confirm that 
both methods let one sample the representative parts along h = and that the STM 
method enables one to obtain the DOS in wider areas. 




Figure 10. Specific heat C/L 2 as a function of T for L = 40. The inset has 
different abscissa and ordinate. 



Because the STM method enables us to obtain the DOS in a wide range of 
sampling space, we can calculate the two-dimensional map of any thermodynamic 
quantity. Figure 12 shows the specific heat and susceptibility per spin as functions of 
T and h when L = 80. This implies that the phase transition temperature converges 
into the Ising case value of 1.1346, as the external field decreases. Related theoretical 
work is found in, e.g., Ref. [42]. 

Figure 13 shows the specific heat as a function of temperature for some values of h 
and L. With positive external field, the phase transition disappears. However, because 
of finite effects, the abnormality, as measured by the diverging behavior, persists to 
some extent. With smaller external field, the divergence behavior remains for larger L. 
Vice versa, with larger L, the more easily it can be shown that the diverging behavior 
disappears. This can be seen as a crossover between L and h. 

Figure 14 shows the dimensionless free energy per spin as a function of 
temperature and external field, which was obtained by applying MBAR to the results 
of the production runs. Note that the partial derivative of this free energy with respect 
to h gives Tfjj, where (M) is defined in (3). The shape at h = suggests a jump of 
m below T c , indicating the existence of first-order phase transitions. 
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Figure 11. Susceptibility x/L 2 as a function of T for L = 40. The inset has 
different abscissa and ordinate. 
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Figure 12. (a) C/L 2 and (b) x/L 2 as functions of T and h for L = 80. The 
solid vertical line corresponds to T = 1.1346, which is the critical temperature of 
the Ising model (in Potts model normalization). 



Finally to study the crossover behavior of the phase transitions, we calculated the 
magnetization by MBAR around the critical point. The scaling form of m is given by 
[43] 

mL^ = ^{tL yt ,hL Vh ), (35) 

where y t = 1/v and = (/3 + According to the crossover scaling formalism 

[43], if t~^th (in Potts model i -14 / 9 /i) is small enough, then the magnetization obeys 
m ~ t 1 / 9 (= t 13 ), and if h~^ht (in Potts model h~ 9 / 1A t) is small enough (i.e., t~ 14 ^ 9 h 
is large enough), then it obeys m ~ h 1 / 14 (= h 1 ^), where t = T y T ° ■ Figure 15(a) 
shows that if the finite-size effects are neglig ible {L 6 / 5 t > 0.1) and t > (l/6)/i 9 / 14 
(i.e., t~ 14 / 9 h is small), then the critical behavior is m ~ t 1 / 9 . Figure 15(b) shows that 
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Figure 13. Specific heat C/L 2 as a function of T for L = 5 (red), 10 (green), 20 
(blue), 40 (magenta), and 80 (cyan), (a) h = 0.0, (b) h = 0.005, (c) h=0.01 (d) 
h = 0.02. 




Figure 14. Free energy per spin f/L 2 and its contour curves as a function of T 
and h. The linear lattice size L was 80. 



if finite-size effects are negligible (L 2S / 15 h > 0.1) and t < (l/6)/i 9 / 14 (i.e., t" 14 / 9 /i is 
large), then the critical behavior is m ~ /i 1 / 14 . Thus, figure 15 clearly shows that the 
line (t = (l/6)/i 9 / 14 ) gives the boundary of the two scaling regimes. 

Because the three-state Potts model in a negative external field is expected to 
behave like the Ising model, we also investigated the crossover behavior between the 
two models. The scaling exponents of Xmecx for increasing L of the Potts model and 
the Ising model are given by j/u = 26/15 and 7/4, respectively. Figure 16 shows 
that the exponents are so similar that we cannot distinguish the difference, despite 
the accuracy of the measurements. 

Thus, we measured different quantities, which are the maximum values of 
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Figure 15. Difference between magnetization and its expected scaling behaviors 
around the critical point. The linear lattice size L was 80. (a) |mL 2 / 15 — 
1.2(L 6 / 5 t) 1 /»| i s illustrated. The solid line corresponds to h = 6i 14 / 9 . (b) 
|mL 2 / 15 - (L 28 / 15 h) 1 / 14 \ is illustrated. The solid line corresponds to h = 6t 14 / 9 . 
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Figures 17-21 show the results. Note that 



and 
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0)/" 



■|max are expected to behave asymptotically as 
respectively, as the lattice size L increases [41]. 
These critical exponents for the Potts model and the two-dimensional Ising model 
are given by v = I and /3 = | so that (1 — = 16/15 and v = 1 and j3 = | 

so that (1 — /3)/f = 7/8 , respectively. We observe that all quantities along h = 
(red curve with filled squares) follow the Potts case, and those with the external field 
at large L (green curve with filled circles and blue curve with filled triangles) follow 
the Ising case. Note that the two curves at h = —0.5 and h = —1.0 converge into 
almost the same line as L increases. On the other hand, the curve at h = —0.5 (green 
curve) is more deviating from the scaling behavior. This can also be understood as 
the crossover between L and h. 
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4. Conclusions 

In this work, we applied Simulated Tempering and Magnetizing (STM) [19, 20] to the 
two-dimensional three-state Potts model. During the simulations, two-dimensional 
random walks in temperature and external field were realized. The random walk 
covered a wide area of temperature and external field so that the STM simulations 
enabled us to study a wide area of the phase diagram from a single simulation run. 

Because of the method's capability of dealing with a wider area of the sampling 
space (as is seen in DOS), we can calculate thermodynamic quantities at an enlarged 
range of the parameter space. We succeeded in reproducing many typical features of 
the system in the presence of an external field. 

We investigated the crossover behavior of the phase transitions by calculating the 
magnetization per spin m around the critical point by reweighting techniques. The 
results showed agreement with previous theoretical studies. Thus, this supports the 
validity of the two-dimensional ST method, or STM. 

With the data of the present work, we can calculate the two-dimensional 
density of states n(E, M) so that we can determine the weight factor for two- 
dimensional multicanonical simulations. Therefore, we can also perform two- 
dimensional multicanonical simulations which would be an interesting future task. 

We also remark that the present methods are useful not only for spin systems but 
also for other complex systems with many degrees of freedom. Note that because this 
method does not require one to change the energy calculations, the method should be 
highly compatible with existing program packages. 
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